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Abstract 



We provide an asymptotically tight, computationally efficient approximation of the joint 
spectral radius of a set of matrices using sum of squares (SOS) programming. The approach 
OO I is based on a search for an SOS polynomial that proves simultaneous contractibility of a finite 

set of matrices. We provide a bound on the quality of the approximation that unifies several 
earlier results and is independent of the number of matrices. Additionally, we present a com- 
I parison between our approximation scheme and earlier techniques, including the use of common 

I quadratic Lyapunov functions and a method based on matrix liftings. Theoretical results and 

• I numerical investigations show that our approach yields tighter approximations. 

1 Introduction 

Stability of discrete linear inclusions has been a topic of major research over the past two decades. 
^ ■ Such systems can be represented as a switched linear system of the form x{k+l) = A„(^i^-jx{k), where 

[ o" is a mapping from the integers to a given set of indices. The above model, and its many variations, 

has been studied extensively across multiple disciplines including control theory, theory of non- 

■ negative matrices and Markov chains, subdivision schemes and wavelet theory, dynamical systems, 
, etc. The fundamental question of interest is to determine whether x{k) converges to a limit, or 

equivalently, whether the infinite matrix products chosen from the set of matrices converge [ BW92t 

■ IDL92t IDLOlj . The research on convergence of infinite products of matrices spans across four 
decades. A majority of results in this area has been provided in the special case of non-negative 
and/or stochastic matrices. A non-exhaustive list of related research providing several necessary 

^ ' and sufficient conditions for convergence of infinite products and their applications includes | CH94t 

IDLOU ILei92l ISWP97j . Despite the wealth of research in this area, finding algorithms that can 
unambiguously decide convergence remains elusive. Much of the difficulty of this problem stems 
from the hardness in computation or efficient approximation of the joint spectral radius of a finite 
set of matrices. This notion was introduced by Rota and Strang [RS60] via the definition 

p{Ai,...,Am) := lira max 1 1 A,,, • • • J | (1) 

fc^oo (Te{l,...,m}'= 

and represents the maximum growth rate that can be achieved by taking arbitrary products of 
the matrices Ai. As in the case of the classical spectral radius, the value of this expression is 
independent of the choice of norm in ([1]). Daubechies and Lagarias |DL92j conjectured that the 
joint spectral radius is equal to a related quantity, the generalized spectral radius, which is defined 
in a similar way except for the fact that the norm of the product is replaced by the spectral radius. 
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Berger and Wang |BW92j proved this conjecture to be true for finite sets of matrices. Blondel and 
Tsitsiklis have shown that computing p is hard from a computational complexity viewpoint, and 
even approximating it is difficult [BTOOal ITB97] . In particular, it follows from their results that 
the problem "Is p < 1?" is undecidable. For rational matrices, the joint spectral radius is not 
a semialgebraic function of the data, thus ruling out a very large class of methods for its exact 
computation. We refer the reader to the survey [BTOObl §3.5] for further results and references on 
the computational complexity of the joint spectral radius. 

It turns out that a necessary and sufficient condition for the stability of a linear difference inclu- 
sion is for the corresponding matrices to have a subunit joint spectral radius, i.e., p{Ai , . . . , Am) < 1; 
see e.g. |SWP97l Thm. 1] and |BT80] . A subunit joint spectral radius is equivalent to the existence 
of a common norm with respect to which all matrices in the set are contractive |Bar88l IKozQOt 
IWir02j : unfortunately, this common norm is in general not finitely constructible. In fact a similar 
result, due to Dayawansa and Martin [DM99j . holds for nonlinear systems that undergo switching. 
A popular approach towards approximating the joint spectral radius or showing that it is indeed 
subunit has been to try to prove simultaneous contractibility (i.e., existence of a common norm 
with respect to which matrices are contractive), by searching for a common ellipsoidal norm, or 
equivalently, a common quadratic Lyapunov function. The benefit of this approach is due to the 
fact that the search for a common ellipsoidal norm can be posed as a semidefinite program and 
solved efficiently using interior point techniques. However, it is not too difficult to generate exam- 
ples where the discrete inclusion is absolutely asymptotically stable, i.e., asymptotically stable for 
all switching sequences, but a common quadratic Lyapunov function, (or equivalently a common 
ellipsoidal norm) does not exist. 

Ando and Shih describe in |AS98j a constructive procedure for generating a set of m matrices 
whose joint spectral radius is equal to but for which no quadratic Lyapunov function exists. 

They prove that the interval [0, is effectively the "optimal" range for the joint spectral radius 
necessary to guarantee simultaneous contractibility under an ellipsoidal norm for a finite collection 
of m matrices. The range is denoted as optimal since it is the largest subset of [0, 1) for which if 
the joint spectral radius is in this subset the collection of matrices is simultaneously contractible 
under an ellipsoidal norm. Furthermore, they show that the optimal joint spectral radius range 
for a bounded set of n x n matrices is the interval [0, -^)- The proof of this fact is based on 
John's ellipsoid theorem [Joh48] . Roughly speaking, John's ellipsoid theorem implies that every 
convex body in n-dimensional Euclidean space that is symmetric with respect to the origin can be 
approximated by inner and outer ellipsoids, up to a factor of Independently, Blondel, Nesterov 
and Theys |BNT05j showed a similar result (also based on John's ellipsoid theorem), that the best 
ellipsoidal norm approximation of the joint spectral radius provides a lower bound and an upper 
bound on the actual value. Given a set M of n x n matrices with joint spectral radius p, and best 
ellipsoidal norm approximation p, it is shown there that 



A major consequence of these results is that finding a common Lyapunov function becomes increas- 
ingly hard as the dimension goes up. 

There have been a number of earlier works proposing different numerical techniques for the 
effective computation of bounds on the joint spectral radius. A natural class of lower bounds is 
obtained by considering periodic switching sequences, in which case only a finite number of matrix 
norms need to be computed. Using a naive approach, the required computational efforts grow 
exponentially as m'^, where k is the period of the sequence. Due to the cyclic property of the 



-^p{M) < p{M) < p{M). 




(2) 
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spectral radius, some terms are redundant, and Maesumi |Mae96| has shown using combinatorial 
techniques that the number of required products can be reduced to /k. Another approach is the 
work of Gripenberg |Gri96j . who has introduced a branch-and-bound algorithm to produce upper 
and lower bounds on the joint spectral radius. Protasov [Pro97t IProOSj has developed a geometric 
method to approximate this quantity, based on a polytopic approximation of a convex set that is 
invariant under the action of the linear operators Ai. This method has also been extended to the 
computation of the so-called p-radius |Pro97] . More recently, Blondel and Nesterov |BN05| have 
proposed an alternative scheme to the computation of the joint spectral radius, by "lifting" the 
matrices using Kronecker products to provide better approximations. A common feature in many 
of these approaches is the presence of convexity-based methods to provide certificates of the desired 
system properties. 

In this paper, we develop a sum of squares (SOS) based scheme for the approximation of the 
joint spectral radius. The method computes, using the techniques of semidefinite programming, a 
homogeneous polynomial that serves as a Lyapunov-like function for the corresponding switched 
linear system. We prove several results on the quality of approximation of the proposed scheme. 
In particular, it will follow from Theorems 13.41 and 14.31 that our SOS-based approximation psos,2d 
satisfies 

r]-2d . psosM^) ^ Pi^) ^ PsosM^)^ 

where rj := min{m, {^~^'^~^)}- To prove this, we use two different techniques, one inspired by recent 
results of Barvinok [Bar02j on approximation of norms by polynomials, and the other one based 
on a convergent iteration similar to that used for Lyapunov inequalities. Our results provide a 
simple and unified derivation of most of the available bounds, including some new ones. We prove 
that the SOS-based approximation is always tighter than that obtained by the use of common 
quadratic Lyapunov functions, and than the one provided by Blondel and Nesterov in |BN05] . 
Furthermore, we show how to compute the bound in [BN05] using matrices that are exponentially 
smaller than those proposed there; this result also follows from the earlier work of Protasov |Pro97] . 
A preliminary version of some of our results has been presented in |PJ07j . 

A description of the paper follows. In Section[2]we present a class of bounds on the joint spectral 
radius based on simultaneous contractivity with respect to a norm, followed by a sum of squares- 
based relaxation, and the corresponding suboptimality properties. In Section [3] we present some 
background material in multilinear algebra, necessary for our developments, and a derivation of a 
bound of the quality of the SOS relaxation. An alternative development is presented in Section SI 
where a different bound on the performance of the SOS relaxation is given in terms of a very natural 
Lyapunov iteration, similar to the classical case. In Section [5] we make a comparison with earlier 
techniques and analyze a numerical example. Finally, in Section [6] we present our conclusions. 

2 Bounds via polynomials and sums of squares 

A natural way of bounding the joint spectral radius is to find a common norm that guarantees 
certain contractiveness properties for all the matrices. In this section, we first revisit this char- 
acterization, and introduce our method of using SOS relaxations to approximate this common 
norm. 

Norms and the joint spectral radius. As we mentioned, there exists an intimate relation- 
ship between the spectral radius and the existence of a vector norm under which all the matrices 
are simultaneously contractive. This is summarized in the following theorem, a special case of 
Proposition 1 in [RS60j by Rota and Strang. 
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Theorem 2.1 ( |RS60| ). Consider a finite set of matrices A = {Ai, . . . , A^}- For any e > 0, there 
exists a norm \\ ■ \\ in M" (denoted as JSR norm hereafter) such that 

ll^ixll < (/j(^) + e) ||x||, Vx G M", i = l,...,m. 

The theorem appears in this form, for instance, in Proposition 4 of |BNT05j . The main idea 
in our approach is to replace the JSR norm that approximates the joint spectral radius with a 
homogeneous SOS polynomial p{x) of degree 2d. As we will see in the next sections, we can 
produce arbitrarily tight SOS approximations, while still being able to prove a bound on the 
resulting estimate. 



Joint spectral radius and polynomials. As the results presented above indicate, the joint 
spectral radius can be characterized by finding a common norm under which all the maps are 
simultaneously contractive. As opposed to the unit ball of a norm, the level sets of a homogeneous 
polynomial are not necessarily convex (see for instance Figured]). Nevertheless, as the following 
theorem suggests, we can still obtain upper bounds on the joint spectral radius by replacing norms 
with homogeneous polynomials. 

Theorem 2.2. Let p{x) he a strictly positive homogeneous polynomial of degree 2d that satisfies 



p{Aix) < 'j'^'^pix) 



, m. 



Then, p{Ai, . . . , Am) < 7- 

Proof. If p{x) is strictly positive, then by compactness of the unit ball in R" and continuity of p{x) 
there exist constants < a < P, such that 
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From the definition of the joint spectral radius in equation ([T]), by taking fcth roots and the limit 
A; — > oo we immediately have the upper bound /o(Ai, . . . , A^) < 7- □ 

The condition in Theorem 12.21 involves positive polynomials, which are computationally hard 
to characterize. A useful scheme, introduced in [ParOOl IPar03] and relatively well-known by now, 
relaxes the nonnegativity constraints to a much more tractable sum of squares (SOS) condition, 
where p{x) is required to have a decomposition as p{x) = J2iPii^)'^- The SOS condition can 
be equivalently expressed in terms of a semidefinite programming (SDP) constraint. In what 
follows, we briefly describe the basic ideas behind SDP and sum of squares programming, and their 
applications to our problem. 
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Semidefinite programming. SDP is a specific kind of convex optimization problem with very 
appealing numerical properties. An SDP problem corresponds to the optimization of a linear 
function over the intersection of an affine subspace and the cone of positive semidefinite matrices. 
For much more information about SDP and its many applications, we refer the reader to the surveys 
|VB961 ITodOl] and the comprehensive treatment in [WSVOOj . 

An SDP problem in standard primal form is usually written as: 

minimize C»X subject to Ai • X = bi, i = 1, . . . ,m 

where C, Ai are symmetric nxn matrices, and X»Y := trace{XY). The symmetric matrix X is the 
optimization variable over which the maximization is performed. The inequality in the second line 
means that the matrix X must be positive semidefinite, i.e., all its eigenvalues should be greater 
than or equal to zero. The set of feasible solutions, i.e., the set of matrices X that satisfy the 
constraints, is always a convex set. In the particular case when C = 0, the problem reduces to 
whether or not the inequality can be satisfied for some matrix X. In this case, the SDP is referred 
to as a feasibility problem. 

There are a number of sophisticated and reliable methods to numerically solve semidefinite 
programming problems. One of the most successful approaches is based on primal- dual interior 
point methods, that generalize many of the techniques used in linear programming [NN94] . The 
interior-point approach to SDP typically involves the iterative solution of a perturbed version of the 
KKT optimality conditions. Each iteration requires the computation of the corresponding Newton 
direction, and the solution of a system of linear equations. A theoretical bound on the number 
of Newton iterations is 0{^/nlog^) for an e-approximate solution. This estimate is signficantly 
more conservative than what is usually experienced in practice, where the dependence on n is 
very mild (typically, 10-40 Newton iterations are enough for most problems). The cost of each 
iteration heavily depends on the structure and sparsity of the matrices Ai, and is dominated by 
the computation of the Hessian and the solution of the corresponding linear system. In the fully 
dense case, this cost is of the order of maxjmn^, m^n^, m^}, where the first two terms correspond 
to the construction of the Hessian, and the last one to the solution of the Newton system. 

Sums of squares programming. Consider a given multivariate polynomial for which we want to 
decide whether a sum of squares decomposition exists. This question is equivalent to a semidefinite 
programming (SDP) problem, because of the following result, that has appeared in different forms in 
the work of Shor j Sho87j . Choi-Lam-Reznick [CLR95j . Nesterov [NesnOj . and Parrilo [ParOOl iParOS] . 

Theorem 2.3. A homogeneous multivariate polynomial p{x) of degree 2d is a sum of squares if 
and only if 

p{x) = {x^'^YQx^'^\ (3) 

where x^'^ is a vector whose entries are (possibly scaled) monomials of degree d in the variables Xi, 
and Q is a symmetric positive semidefinite matrix. 

Since in general the entries of j;''^! are not algebraically independent, the matrix Q in the 
representation ([3|) is not unique. In fact, there is an affine subspace of matrices Q that satisfy 
the equality, as can be easily seen by expanding the right-hand side and equating term by term. 
To obtain an SOS representation, we need to find a positive semidefinite matrix in this affine 
subspace. Therefore, the problem of checking if a polynomial can be decomposed as a sum of 
squares is equivalent to verifying whether a certain affine matrix subspace intersects the cone of 
positive definite matrices, and hence an SDP feasibility problem. 
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Example 2.4. Consider the quartic homogeneous polynomial in two variables described below, and 
define the vector of monomials as [x"^ , y"^ , xyf'" . 

p{x, y) = 2x^ + 2x^y — x^y^ + 5y^ 



x^ 
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. ^y . 



= Q'lix^ + 922/ + (933 + 2912)2:^1/^ + 2qi3X^y + 2q23xy^ 

For the left- and right-hand sides to be identical, the following linear equations should hold: 

911=2, ^22 = 5, g33 + 2gi2 = -l, 2gi3 = 2, 2^23 = 0. (4) 

A positive semidefinite Q that satisfies the linear equalities can then be found using SDP. A 
particular solution is given by: 
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and therefore we have the sum of squares decomposition: 



p{x,y) = ^(2x2 



3y^ + xyf + -{y^ + 3xyf. 



□ 



PSOS,2d : = 



inf 7 

p(2:)6lR2dH,7 



S.t. 



2.1 Norms and SOS polynomials 

The procedure described in the previous subsection can be easily adapted to the case where the 
polynomial p{x) is not fixed, but instead we search for an SOS polynomial in a given affine family 
(for instance, all homogeneous polynomials of a given degree). 

This line of thought immediately suggests the following SOS relaxation of the conditions in 
Theorem O 

p{x) is SOS , V 

^2dp^x) - p{Aix) is SOS ^ ' 

where K2d[x] is the set of homogeneous polynomials of degree 2d. 

Remark 2.5. Theorem \2.2\ requires a strictly positive polynomial p{x), so it would be natural to 
add some strict positivity condition to the relaxation ^B^. For instance, one could require for the 
polynomial p{x) to belong to the relative interior of the SOS cone. However, since interior-point 
methods by construction always produce solutions in the relative interior of the corresponding convex 
set, this is automatically satisfied if the problem is feasible. Alternatively, it is possible to give a 
formulation that includes terms of the form e\\x\\'^'^, for small positive e. These modifications are 
unnecessary in practice. 

For any fixed degree d and any given 7, the constraints in this problem are all of SOS type, 
and thus equivalent to semidefinite programming. Therefore, the computation of psos,2d is a 
quasiconvex problem, and can be easily solved with a standard SDP solver, and a simple bisection 
method for the scalar variable 7. By Theorem 12. 2 1 the solution of this relaxation yields an upper 
bound on the joint spectral radius 



p{Ai, . . .,Am) < PSOS,2d, 

where 2d is the degree of the approximating polynomial. 



(6) 
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2.2 Quality of approximation 

What can be said about the quality of the bounds produced by the SOS relaxation? We present next 
some results to answer this question; a more complete characterization is developed in Section 13.11 
An inspiring result in this direction is the following theorem of Barvinok, that quantifies how tightly 
SOS polynomials can approximate norms: 

Theorem 2.6 ( |Bar02| . p. 221). Let \ \ ■ \ \ be a norm in W\ For any integer d > 1 there exists a 
homogeneous polynomial p{x) in n variables of degree 2d such that 

1. The polynomial p{x) is a sum of squares. 

2. For all x E M", 

p{x)2d < ||x|| < k{n,d) p{x)2d ^ 

where k{n,d) := ("+;^~^)^. 

For fixed state dimension n, by increasing the degree d of the approximating polynomials, the 
factor in the upper bound can be made arbitrarily close to one. In fact, for large d, we have the 
approximation 

To apply these results to our problem, consider the following. If p{Ai, . . . ,Am) < 7, by Theo- 
rem [27T] (and sharper results in |Bar881 IKozQOt [Wir02| ) there exists a norm || • || such that 

I I < 7I l^l I, Vx G M", i = 1, . . . , m. 

By Theorem 12.61 we can therefore approximate this norm with a homogeneous SOS polynomial 
p{x) of degree 2d that will then satisfy 

p{Aix)^ < \ \Aix\\ < 7||3;|| < ^ k{n,d) p{x)'2d ^ 

and thus we know that there exists a feasible solution of 

J p{x) is SOS 

\ a'^'^p{x) — p{Aix) >0 i = l,...,m, 

for a = k{n, d)p{Ai, . . . , Am)- 

Despite these appealing results, notice that in general we cannot yet conclude from this that 
the proposed SOS relaxation will always obtain a solution that is within k{n,d)~^ from the true 
spectral radius. The reason is that even though we can prove the existence of a p{x) that is SOS 
and for which a'^'^p{x) —p{Aix) are nonnegative for all i, it is unclear whether the last m expressions 
are actually SOS. We will show later in the paper that this is indeed the case. Before doing this, 
we concentrate first on two important cases of interest, where the described approach guarantees a 
good quality of approximation. 

Planar systems. The first case corresponds to two-dimensional (planar) systems, i.e., when 
n = 2. In this case, it always holds that nonnegative homogeneous bivariate polynomials are SOS 
(e.g., [RezOO] ). Thus, we have the following result: 

Theorem 2.7. Let {Ai, . . . , Am} C M^^^. Then, the SOS relaxation always produces a solution 
satisfying: 

\PSOS,2d <{d+ 1)~ 2d PSOS,2d < P{Al, Am) < PSOS,2d- 

This result is independent of the number m of matrices. 
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Figure 1: Level sets of the quartic homogeneous polynomial V{xi,X2)- These define a Lyapunov 
function, under which both Ai and A2 are (1 + e)-contractive. The value of e is here equal to 0.01. 

Quadratic Lyapunov functions. In the quadratic case (i.e., 2d = 2), it is also true that 
nonnegative quadratic forms are sums of squares. Since 



d VI 



n, 



the inequality 

PSOS,2 < p{Ai, . . . , Am) < PSOS,2 (7) 



follows. This bound exactly coincides with the results of Ando and Shih |AS98j or Blondel, Nesterov 
and Theys |BNT05| . This is perhaps not surprising, since in this case both Ando and Shih's 
proof [AS98] and Barvinok's theorem rely on the use of John's ellipsoid to approximate the same 
underlying convex set. 

Level sets and convexity Unlike the norms that appear in Theorem 12.11 an appealing feature 
of the SOS-based method is that we are not constrained to use polynomials with convex level sets. 
This enables in some cases much better bounds than what is promised by the theorems above, as 
illustrated in the following example. 

Example 2.8. This is based on a construction by Ando and Shih IAS98^ . Consider the problem of 
proving a bound on the joint spectral radius of the following matrices: 



A, 



1 
1 



Ao 



1 

-1 



For these matrices, it can be easily shown that p{Ai, A2) = 1. Using a common quadratic Lyapunov 
function (i.e., the case d = 2), the upper bound on the joint spectral radius is equal to \f2. However, 
a simple quartic SOS Lyapunov function is enough to prove an upper bound ofl + e for every e > 0, 
since the SOS polynomial 

V{x) = {xj - xlf + e{xl + xlf 
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satisfies 



il + e)Vix)-V{Aix) = ixl-xl + e{xl + xl)f 
{l + e)V{x)-ViA2x) = {xl-xl + eixl + xDf. 

The corresponding level sets ofV{x) are plotted in FigureUl o.'^^d are clearly non-convex. 

3 Symmetric algebra and induced matrices 

We present next some further bounds on the quahty of the SOS relaxation ([5]), either by a more 
refined analysis of the SOS polynomials in Barvinok's theorem or by explicitly producing an SOS 
Lyapunov function of guaranteed suboptimality properties. These constructions are quite natural, 
and parallel some lifting ideas as well as the classical iteration used in the solution of discrete-time 
Lyapunov inequalities. Before proceeding further, we briefly revisit some classical notions from 
multilinear algebra. 

Symmetric algebra of a vector space Consider a vector x G M", and an integer d > 1. We 
define its d-lift x^'^'^ as a vector in M^, where N := ("'^^~^), with components {y/a^. x°'}a, where 



a = (ai, . . . , a„), \a\ := ai = d, and a! denotes the multinomial coefficient a! := 



J 



i. That is, the components of the lifted vector are the monomials of degree d, scaled by 



the square root of the corresponding multinomial coefficients. 
Example 3.1. Let n = 2, and x = [u,v]'^ . Then, we have 
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.,3 



The main motivation for this specific scaling of the components, is to ensure that the lifting 
preserves some of the properties of the underlying normed space. In particular, if 1 1 • 1 1 denotes the 
standard Euclidean norm, it can be easily verified that = Thus, the lifting operation 

provides a norm-preserving (up to power) embedding of M" into M^. When the original space is 
projective, this is the so-called Veronese embedding. 

This concept can be directly extended from vectors to linear transformations. Consider a linear 
map in M", and the associated nx n matrix A. Then, the lifting described above naturally induces 
an associated map in M^, that makes the corresponding diagram commute. The matrix representing 
this linear transformation is the d-th induced matrix of A, denoted by A^"^, which is the unique 
N X N matrix that satisfies 



[d] 



In systems and control, these classical constructions of multilinear algebra have been used under 
different names in several works, among them |Bro741 [Zel94] and (imphcitly) |BN05] . Although not 
mentioned in the Control literature, there exists a simple explicit formula for the entries of these 
induced matrices; see |Mar73t IMM92| . The d-th induced matrix A^'^ has dimensions N x N. Its 
entries are given by 

{Al^^U = ^^^tl^, (8) 
V^/i(a)^(/3) 
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where the indices a, (3 are all the d-element multisets of {1, . . . ,n}, the notation per indicates the 
permanent of a square matrix, and fJ-{S) is the product of the factorials of the multiplicities of the 
elements of the multiset S. 

Example 3.2. Consider the case n = 2, d = 3. The corresponding 3-element multisets are {1, 1, 1}, 
{1,1,2}, {1,2,2} and {2,2,2}. The third induced matrix is then 

afi V3aiiai2 \/3aiiai2 "12 

^[3]^ \/3afia2i aii(aiia22 + 2021012) 012(2011022 + 021012) \/3of2022 

a/30ii02i 021(2011022+021012) 022(011022 + 2021012) \/3oi2022 
O21 \/302i022 \/302l022 "22 

It can be shown that these operations define an algebra homomorphism, i.e., they respect the 
structure of matrix multiplication. In particular, for any matrices A, B of compatible dimensions, 
the following identities hold: 

{AB)'^'^^ = A^<^^B^'^\ (A-^)['^l = (^['^1)-^ 

Furthermore, there is a simple and appealing relationship between the eigenvalues of A^'^l and those 
of A. Concretely, if Ai,...,An are the eigenvalues of A, then the eigenvalues of ^'^^ are given 
by njgs"^j where S C {1, . . . , n}, jS"! = d; there are exactly ("''^^~^) such multisets. A similar 
relationship holds for the corresponding eigenvectors. Essentially, as explained below in more 
detail, the induced matrices are the symmetry-reduced version of the d-fold Kronecker product. 

The symmetric algebra and associated induced matrices are classical objects of multilinear alge- 
bra. Induced matrices, as defined above, as well as the more usual compound matrices, correspond 
to two specific isotypic components of the decomposition of the d-fold tensor product under the 
action of the symmetric group S"^ (i-e., the symmetric and skew- symmetric algebras). Compound 
matrices are associated with the alternating character (hence their relationship with determinants) , 
while induced matrices correspond instead to the trivial character, thus the connection with per- 
manents. Similar constructions can be given for any other character of the symmetric group, by 
replacing the permanent in ([8]) with the suitable immanants; see |Mar73| for additional details. 

3.1 Bounds on the quality of psos,2d 

In this section we present a bound on the approximation properties of the SOS approximation, 
based on the ideas introduced above. As we will see, the techniques based on the lifting described 
will exactly yield the factor k{n,d)~^ suggested by Barvinok's theorem. 

We first prove a preliminary result on the behavior of the joint spectral radius under d-lifting. 
The scaling properties described earlier can be applied to obtain the following: 

Lemma 3.3. Given matrices {Ai, . . . ,Am} C M"^" and an integer d > 1, the following identity 
holds: 

p(4"],...,aM) = p(^i,...,A^)'^. 

The proof follows directly from the definition ([T]) and the two properties {AB)^'^ = A^'^B^'^, 
\\x^'^\ \ = and it is thus omitted. 

Combining all these inequalities, we obtain the main result of this paper: 

^The permanent of a matrix A G R"'^" is defined as per(yl) := X^^gn YYi=i where n„ is the set of all 

permutations in n elements. 
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Theorem 3.4. The SOS relaxation satisfies: 



("^d"') " Psos,2d < p{Ai, ...,Am)< PsosM- (9) 



Proof. Since the dimension of is C^^^ from Lemma 13.31 and inequahty ([7]) it fohows that: 

Combining this with ([6]) and the inequahty (proven later in Theorem I5.ip , 

PsosMAu • • • , Amf < PsosM?.. . . , ^[^1), 
the result follows. □ 



4 Sum of squares Lyapunov iteration 

We describe next an alternative approach to obtain bounds on the quality of the SOS approximation. 
As opposed to the results in the previous section, the bounds now explicitly depend on the number 
of matrices, but will usually be tighter in the case of small m. 
Consider the iteration defined by 

^ m 

Vo{x) = 0, Vk+i{x) = Q{x) + - ^ VkiAix), (10) 

^ i=i 

where Q{x) is a fixed n-variate homogeneous polynomial of degree 2d and /3 > 0. The iteration 
defines an affine map in the space of homogeneous polynomials of degree 2d. As usual, the iteration 
will converge under certain assumptions on the spectral radius of this linear operator. 

Theorem 4.1. The iteration defined in U0\) converges for arbitrary Q{x) z//9(Af'^^+- • •+^^'^') < [3. 

Proof. The vector space of homogenous polynomials M2d[xi, . . . , x^] is naturally isomorphic to 
the space of linear functionals on (M")[2'^l, via the identification Vk{x) = (w^, xt^'^l), where G 

m( 2d ) is the vector of (scaled) coefficients of Vk{x). Then, since Vk{Aix) = (u^, (j4j2;)[^'^l) = 
(ufc, 74p'^'x[^'^l) = ((A|^'^^)-^Ufc, the iteration pUj) can be simply expressed as: 



T 

Vk, 



and it is well known that an affine iteration converges if the spectral radius of the linear term is 
less than one. □ 

For simplicity of notation, we define the following quantity, corresponding to the spectral radius 
of the sum of the 2ci- lifted matrices: 

PSR,2d:=pi4'^+--- + Ag^)i-.. (11) 
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Theorem 4.2. The following inequality holds: 

PSOS,2d < PSR,2d 

Proof. Choose a Q{x) that is in the interior of the SOS cone, e.g., Q{x) := {Yli=i ^fY-, l^t 
/3 = p{Al^'^'^ -\- ■ ■ ■+ALm'^)+e. The iteration guarantees that V^+i is SOS if is. By induction, ah 
the iterates Vk are SOS. By the choice of (5 and Theorem l4.lt the Vk converge to some homogeneous 
polynomial Voo{x). By the closedness of the cone of SOS polynomials, the limit Voo is also SOS. 
Furthermore, we have 

PV^{x) - V^{Aix) = l3Q{x) + ^ Voo{A,x) 

and therefore the expression on the left-hand side is SOS. This implies that p{x) := V^(x) is a 
feasible solution of the SOS relaxation (0). Taking e ^ 0, the result follows. □ 

Notice that if the spectral radius condition in Theorem 14.11 is satisfied, then for any fixed Q{x) 
the corresponding limit T4o(x) = {voo,x^'^'^) can be simply obtained by solving the nonsingular 
system of linear equations 

thus generalizing the standard Lyapunov equation. The iteration argument is only used to prove 
that the solution of this linear system yields a strictly positive SOS polynomial. A slightly different 
approach here is via the finite-dimensional version of the Krein-Rutman theorem (or generalized 
Perron- Probenius); see for instance |Pro97j or [PKOOj . 

Theorem 4.3. The SOS relaxation satisfies: 

m'2d psos,2d < • • • , Am) < Psos,2d- 

Proof. This follows directly from inequality ([6]), and the fact that 

< ■ PiA?'^ , . . . , ^P'^l)^ = • p {A,, . . . , A^) , 

where the first inequality is Theorem 14.21 the second one follows from the general fact that p{Ai + 
• • • + Am) < mp{Ai, . . . , Am) (see e.g.. Corollary 1 in [BN05] ). and the third from Lemma [331 □ 

The iteration (jlOh is the natural generalization of the Lyapunov recursion for the single matrix 
case, and of the construction by Ando and Shih in |AS98j for the quadratic case. By the remarks 
in Section [3] above, and as described in more detail in the next section, it can be shown that the 
quantity psR,2d is essentially equal to those defined by Protasov in [Pro971 §4] and Blondel and 
Nesterov in |BN05| . As a consequence of Theorem l4.2l the SOS-based approach will always produce 
estimates at least as good as the ones given by these procedures. 

5 Comparison with earlier techniques 

In this section we compare the psos,2d approach with some earlier bounds from the literature. We 
show that our bound is never weaker than those obtained by all the other procedures. 
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Steps / 2d 


Accuracy 


[BN05J, Kronecker 


[BN05J 


, semidefinite 


This paper 


n = 2 


n = 10 


n = 2 


n = 10 


n = 2 


n = 10 


1/2 


0.707 


4 


100 


3 


55 


3 


55 


2/4 


0.840 


16 


10000 


6 


1540 


5 


715 


3/8 


0.917 


256 


108 


21 


1186570 


9 


24310 


4/16 


0.957 


65536 


1016 


231 


7.04 X IQii 


17 


2042975 


5/32 


0.978 


4.29 X 10^ 




26796 


2.48 X 10^3 


33 


3.5 X 10^ 



Table 1: Comparison of matrix sizes for the different Ufting procedures to compute psR,2d- The 
matrix size for the Kronecker hfting is n?'^, while the recursive semidefinite lifting is given by the 
d-step recursion S2k = ('^''2'^) with si = n, and the size for the symmetric algebra approach is 
(""'"m"^)' '^^^ accuracy estimates correspond to the case of two matrices, i.e., m = 2. 



5.1 Methods of Protasov and Blondel-Nesterov 

Protasov |Pro97j has shown that an upper bound on the "standard" joint spectral radius can be 
computed via the so-called joint p-radius, a generalization of the definition ([1]) involving p-norms. 
Furthermore, he has shown that in the case of even integer p, the value of the p-radius of an 
irreducible finite set of matrices exactly corresponds to the spectral radius of a single operator, 
that can in principle be constructed based on the matrices Ai. 

Independently, Blondel and Nesterov |BN05] developed a technique based on the calculation 
of the spectral radius of "lifted" matrices. In fact, they present two different lifting procedures 
("Kronecker" and "semidefinite" liftings), and in Section 5 of their paper, they describe a family 
of bounds obtained by arbitrary combinations of these two liftings. 

Both of these methods are in fact equivalent to our construction of psR,2d in Section [U in 
the sense that they all yield exactly the same numerical value. By Theorem 14. 2| they are thus 
also weaker than the SOS-based construction. The bound defined by psR,2d in (fTT|l relies on a 
single canonically defined lifting, and requires much less numerical effort than the Blondel-Nesterov 
construction. Furthermore, instead of the somewhat more complicated construction of Protasov, 
the expression of the entries of the lifted matrices are given by the simple formula dS]), making a 
computer implementation straightforward, with no irreducibility assumptions being required. 

It can be shown that our construction (or Protasov's) exactly corresponds to a fully symmetry- 
reduced version of the Blondel-Nesterov procedure, thus yielding equivalent bounds, but at a much 
smaller computational cost since the corresponding matrices are exponentially smaller (for fixed n, 
the size grows as 0(d"~i) as opposed to 0{n'^'^)). Therefore, even if no SDPs are to be solved (as 
would be required by the tighter bound psos,2d), the formulation in terms of the matrices A^^^ 
still has many advantages. 

As an illustrative comparison of the advantages of this reduced formulation, in Table [1] we 
present the sizes of the matrices required by the method in |BN05j (using the "Kronecker" and 
"recursive semidefinite" liftings) and our approach to psR,2d via the symmetric algebra. The data 
in Table [T] corresponds to that in |BN051 p. 266] (with a minor misprint corrected). 
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5.2 Common quadratic Lyapunov functions 

This method corresponds to finding a common quadratic Lyapunov function, either directly for the 
matrices Ai, or for the lifted matrices A^f^. Specifically, let 

PCQ,2, ■■= inf { 7 I l''P - {A^YpA? to, Pyo}. 

This is essentially equivalent to what is discussed in Corollary 3 of |BN05j , except that the matrices 
involved in our approach are exponentially smaller (of size ("^^~^) rather than n'^), as all the 
symmetries have been taken ouo Notice also that, as a consequence of their definitions, we have 

We can then collect most of these results in a single theorem: 
Theorem 5.1. The following inequalities between all the bounds hold: 

p{Ai, Am) < PSOS,2d < PCQ,2d < PSR,2d- (12) 



Proof. The left-most inequality is Q. The right-most inequality follows from a similar (but 
stronger) argument to the one given in Theorem 14.21 above, since the spectral radius condition 
p[Af^ _j_ . . . _j_ j4[^'^1) < p actually implies the convergence of the matrix iteration in 5^ given by 

i=i 

For the middle inequality, let p{x) := {x^'^)'^ Px^'^ . Since P ;^ 0, it follows that p{x) is SOS. 
From — (A^f^)'^ PA^^ >z 0, left- and right-multiplying by x^'^'^, we have that j'^'^p{x) — p{Aix) 
is also SOS, and thus p{x) is a feasible solution of ([5|), from where the result directly follows. □ 

Remark 5.2. We always have psos,2 = PCQ,2, since both correspond to the case of a common 
quadratic Lyapunov function for the matrices Ai . 

5.3 Computational cost 

In this section we quantify the computational cost of the bound psos,2d- In the following calculations 
we keep d fixed, and study the scaling behavior as a function of the dimension n. 

As mentioned in Section [21 solving a semidefinite programming problem typically requires sev- 
eral Newton iterations, with the cost of each iteration being dominated by the construction of the 
Hessian and solution of the corresponding linear system. For the SOS bound psos,2d, the underly- 
ing SDP problem has m + 1 matrix inequalities corresponding to the SOS constraints in each 
of dimension (^~^'^~^) ~ ^ • n'^, which is 0{n'^) for fixed d. The number of decision variables is 
approximately m ■ (""''2^"^) ^ " ^^'^^ Thus, using a simple bisection method for 7, exploiting 
the block-diagonal structure, and the fact that the number of Newton iterations is essentially con- 
stant, we obtain that the approximate cost of obtaining an e-approximate solution of psos,2d is 
0{m ■ n^'^ ■ log i), where d is chosen such that e ~ f^^x^ or e ~ m~2d^ depending on whether we 
use bounds that depend on the number of matrices (Theorem 14. 3p or not (Theorem 13. 4p . 

^There seems to be a typo in equation (7.4) of [BN05] . as all the terms should likely read Af". 
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d 






PSOS,2d 


PCQ,2d 


PSR,2d 


1 


4 


10 


9.761 


9.761 


12.519 


2 


10 


35 


8.92 


9.01 


9.887 


3 


20 


84 


8.92 


8.92 


9.3133 



Table 2: Comparison of the different approximations for Example 15.41 



We remark that these quantities are a relatively coarse estimate of the best possible algorithmic 
complexity, since very little structure of the corresponding SDP problem is being exploited. It 
is known that for structured problems such as the ones appearing here much more efficient SDP- 
based algorithms can be developed. In particular, in the context of sum of squares problems several 
techniques are known to exploit some of the available structure for more efficient computation; see 
|(mND03[ [LPM IRV06] . 

5.4 Examples 

We present next two numerical examples that compare the described techniques. In particular, we 
show that the bounds in Theorem 15.11 can all be strict. 

Example 5.3. Here we revisit the construction presented earlier in Example \2.8[ For the matrices 
given there we have: 



PSOS,2 
PSOSA 



V2, 

1, 



PCQ,2 
PCQA 



1. 



PSR,2d 



Example 5.4. Consider the three 4x4 matrices (randomly generated) given by: 



A, 



1 

6 
-1 




7 
-2 
-2 

9 



4 
-3 
-6 

1 



-3 
-2 
4 
1 



3 
1 

-3 
-5 




4 
1 

-1 



A. 



1 


-1 



4 
5 
-1 
5 



10 
-4 
6 
1 



The value of the different approximations are presented in TablelB A lower bound is ^(^1^3) 2 
8.9149, which is extremely close (and perhaps exactly equal) to the upper bound psoSA- Notice from 
the d = 2 entry of Table [H that all the inequalities flgj) can be strict. 



6 Conclusions 

We introduced a novel scheme for the approximation of the joint spectral radius of a set of matrices 
using sum of squares programming. The method is based on the use of a multivariate polynomial to 
provide a norm-like quantity under which all matrices are contractive. We provided an asymptoti- 
cally tight estimate for the quality of the bound, which is independent of the number of matrices. 
We also proposed an alternative bound, that depends on the number m of matrices, based on a 
generalization of a Lyapunov iteration. 

Our results can be alternatively interpreted in a simpler way as providing a trajectory-preserving 
lifting to a higher dimensional space, and proving contractiveness with respect to an ellipsoidal 
norm in that space. In this case, a weaker estimate can be obtained by computing the spectral 
radius of a fixed matrix. These results generalize earlier work of Ando and Shih [AS98] . Blondel, 
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Nesterov and Theys |BNT05j . and provide an improvement over the lifting procedure of Blondel 
and Nesterov [BN05| . The good performance of our procedure was also verified using numerical 
examples. 
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